Saumya Biswas (saumyab@uoregon.edu)
In a periodic crystal with periodic boundary condition, the wave function of a particle assummes the bloch form resolvable into the product of a plain wave($|k \rangle$) and a cell periodic part($| u_{n}(k) \rangle$). \begin{eqnarray} |\psi_n(k) \rangle = |k \rangle \otimes | u_{n}(k) \rangle \end{eqnarray}
The periodicity of the crystal limits the crystal momenta, k to a set of discrete values(as many as the number of cells(L) in the crystal). We name this set $k_{lat}$ and enumerate the crystal momenta of the $k_{lat}$ with the wave-vector associated with them.
In a 1d crystal of L cells, for $k \in k_{lat}$, k is one of the values, \begin{eqnarray} k = \frac{2\pi}{La}\times n, \ \ \ \ \ \text{where} \ \ \ n = 0,1,2,...., L-1 \ \ \text{and a is the length of the unit cell} \end{eqnarray}
We enumerate the position operator with the index, j. So a complete set of position eigenkets have corresponding position eigenvalues labeled as, \begin{eqnarray} x = ja ,\ \ \ \ \ \ \ \text{where} j = 0, 1, 2, ..., L-1 \end{eqnarray} Each position eigenket defined as such simply refer to a particle localized at a certain cell, j. We label these position eigenkets as $|j\rangle$
The momentum eigenfunctions in the position basis are the properly normalized plane waves(continuum state normalization). \begin{eqnarray} \langle j|k \rangle = \sqrt{\frac{a}{h}} e^{i \frac{2\pi x}{L}j } \\ \langle k'| k \rangle = \delta(k-k') \end{eqnarray} The $\langle j|k \rangle$ have the dimension of $[\text{momentum dimension}]^{-\frac{1}{2}}$
The matrix elements of the crystal momentum operator($\hat{\bf{k}}$) in the basis of ${|j \rangle}$ are formed as follows \begin{eqnarray} \langle j|k|j' \rangle = \sum\limits_{k \in k_{lat}} \langle j | k \rangle \langle k |j' \rangle \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\frac{1}{h} \frac{2\pi}{L}\sum\limits_{n=0}^{L-1} n e^{-i \frac{2\pi n}{L}(j'-j) } \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\frac{2\pi}{h}\left\{ \begin{array}{ll} (L-1)/2, \ \ \text{when} \ \ \ j=j'\\ \frac{1}{e^{i\frac{2\pi}{L}(j-j')}-1}, \ \text{otherwise} \end{array} \right. \end{eqnarray}
In [48]:
from qutip import *
import matplotlib.pyplot as plt
import numpy as np
In [49]:
periodic_atom_chain8 = Lattice1d(num_cell=8, boundary = "periodic")
k8 = periodic_atom_chain8.k()
In [50]:
[ks8, pw8] = k8.eigenstates()
In [51]:
ks8 # In units of 2*pi/(L*a), if ks[1] = 1, the wavevector/crystal-momentum of the
# corresponding eigen-vector of k is a plane wave with wavelength of L*a, the length
# of the crystal. ks[2] = 2 has a corresponding plane wave eigenvector with
# wavelength of 2*L*a, twice the length of the crystal.
Out[51]:
In [70]:
pw8[0] # ks[0] = 0 has an eigenvector that has a frequency of 0 and wavelength of
# infinity, the eigen-vector is normalized to unity, meaning the sum of all
# the elements sqared is 1
Out[70]:
In [53]:
# The eigenvectors are better demonstrated with a latice of longer length
# Here we plot the second and third eigen-vectors(real and imaginary parts) of the first
# 3 eigenvectors of the crystal momentum operator for a 1d lattice with 64 cells
periodic_atom_chain64 = Lattice1d(num_cell=64, boundary = "periodic")
k64 = periodic_atom_chain64.k()
[ks64, pw64] = k64.eigenstates()
plt.plot(np.real(pw64[1]),'r')
plt.plot(np.imag(pw64[1]),'b')
plt.plot(np.abs(pw64[1]),'g')
plt.xlabel('position')
plt.ylabel('crystal momentum eigen vector')
plt.show()
plt.close()
In [54]:
k64 = periodic_atom_chain64.k()
[ks64, pw64] = k64.eigenstates()
plt.plot(np.real(pw64[2]),'r')
plt.plot(np.imag(pw64[2]),'b')
plt.plot(np.abs(pw64[2]),'g')
plt.xlabel('position')
plt.ylabel('crystal momentum eigen vector')
plt.show()
plt.close()
In [55]:
Hamt8 = periodic_atom_chain8.Hamiltonian()
In [56]:
pw8_M = np.array([pw8[0].full(),pw8[1].full(),pw8[2].full(),pw8[3].full(),
pw8[4].full(),pw8[5].full(),pw8[6].full(),pw8[7].full()])
pw8_M = np.squeeze(pw8_M, axis=2) # changing shape from (8,8,1) to (8,8)
# pw8_M is a matrix of 8 columns, each of which are the eigenvectors
# of the crystal momentum operator
pw8_M = Qobj(pw8_M)
In [57]:
pw8_M * Hamt8 * pw8_M.dag()
Out[57]:
So, eigenvectors of the crystal momentum are eigenvectors of the Hamiltonian as well. This is a consequence of the translational symmetry of the lattice. The generator of the lattice translational operator is the crystal momentum operator. Due to the translational symmetry, the lattice translational operator, the crystal momentum and the Hamiltonian all commute with each other. So they can have simultaneous eigenstates.
In [58]:
# A check that the Hamiltonian and the crystal momentum operator do indeed commute
Hamt8 * k8 - k8 * Hamt8
Out[58]:
We form a Gaussian state in the position basis and evaluate its Discrete Fourier Transform(DFT) for evaluating its frequency components.
In [112]:
xs8 = np.arange(0,8)
Gaussian8_sum_to_1 = 1/np.sqrt(2*np.pi* 0.5 **2) * np.exp(-(xs8 - 3)**2/2/0.5**2)
Gaussian_state8 = np.sqrt(Gaussian8_sum_to_1) # A normalized Gaussian state
In [113]:
plt.plot(Gaussian_state8,'r')
plt.xlabel('position')
plt.ylabel('Probability Amplitude')
plt.show()
plt.close()
In [115]:
k_exps = pw8_M.dag() * Gaussian_state # ket vector of inner products with crystal
# momentum eigenkets
k_exps = k_exps/ np.sqrt( np.sum(np.multiply(k_exps, np.conj(k_exps)))) # normalization
In [116]:
dft_Gaussian_state = np.fft.fft(Gaussian_state)
shf_dft_Gaussian_state = np.roll(dft_Gaussian_state,4)
sq_dft_Gaussian_state = np.multiply(shf_dft_Gaussian_state, np.conj(shf_dft_Gaussian_state))
norml_dft_Gaussian_state = shf_dft_Gaussian_state / np.sqrt(np.sum(sq_dft_Gaussian_state))
In [118]:
plt.plot(np.abs(norml_dft_Gaussian_state),'b')
plt.plot(np.abs(k_exps),'ro')
plt.xlabel('frequency number')
plt.ylabel('frequency Amplitude')
plt.show()
plt.close()
Therefore, expanding a state in momentum eigenkets is equivalent of a DFT for a periodic lattice.
In [119]:
qutip.about()
In [120]:
qutip.cite()